Hydrogen Dynamics in Superprotonic CSHSO4 
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We present a detailed study of proton dynamics in the hydrogen-bonded superprotonic conductor 
CsHSCU from first-principles molecular dynamics simulations, isolating the subtle interplay between 
the dynamics of the O-H chemical bonds, the O- ■ ■ H hydrogen bonds, and the SO4 tetrahedra in 
promoting proton diffusion. We find that the Grotthus mechanism of proton transport is primarily 
responsible for the dynamics of the chemical bonds, whereas the reorganization of the hydrogen-bond 
network is dominated by rapid angular hops in concert with small reorientations of the SO4 tetra- 
hedra. Frequent proton jumping across the O-H- • • O complex is countered by a high rate of jump 
reversal, which we show is connected to the dynamics of the SO4 tetrahedra, resulting in a diminished 
CsHSCU/CsDSCU isotope effect. We also find evidence of multiple timescales for SO4 reorientation 
events, leading to distinct diffusion mechanisms along the different crystal lattice directions. Finally, 
we employ graph-theoretic techniques to characterize the topology of the hydrogen-bond network 
and demonstrate a clear relationship between certain connectivity configurations and the likelihood 
for diffusive jump events. 

PACS numbers: 71.15.Pd,66.30.-h,66.30.Dn,82.47.Pm 
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I. INTRODUCTION AND MOTIVATION 

Despite extensive efforts to realize the hydrogen econ- 
omy, key technological hurdles to the widespread adop- 
tion of fuel-cell technology remain. One of the sev- 
eral challenges is the discovery and optimization of 
proton-conducting electrolyte materials for mid- to high- 
temperature use. Recently, a great deal of research ef- 
fort has focused on anhydrous solid-state materials, in 
part because these offer greater flexibility in their range 
of operating temperatures than materials containing liq- 
uid water. Among the promising candidates are vari- 
ous derivatives of CSHSO4, which has been shown to 
exhibit high ionic conductivity (> 10~ 2 (Sl-cm)" 1 )^ at 
target operating temperatures. Fuel-cell operation us- 
ing electrolytes based upon CSHSO4 and similar mate- 
rials has already been successfully demonstrated in the 
laboratory^ 3 - A complete theoretical picture of the de- 
tailed atomistic mechanisms involved in proton transport 
in these materials is highly desirable, since such knowl- 
edge would be useful in directing future research efforts in 
the optimization and adaptation of new materials. Previ- 
ous molecular dynamics studies^iS of CSHSO4 based on 
fitted interatomic potentials have aided in highlighting 
the basic phenomenology of proton transport, but such 
investigations are unable to capture the full complexity 
of hydrogen bonding and electronic interactions, particu- 
larly in a dynamic environment that features rapid bond 
breaking and forming. The work of Ke and Tanaka 7 -^ in- 
corporated first-principles methodologies, but their anal- 
ysis was grounded in static rather than dynamics cal- 
culations. The present study aims to elucidate the de- 
tailed atomistic pathways and mechanisms involved in 
hydrogen diffusion in superprotonic CSHSO4 using first- 
principles molecular dynamics. 

CSHSO4 was the first known crystalline material to 
exhibit both hydrogen bonding and superprotonic be- 



havior and has among the higher ionic conductivities of 
the known solid- acid materials. The room-temperature 
phase is monoclinic and features a static, well-defined 
network of hydrogen bonds. The higher-temperature su- 
perprotonic phase (usually designated Phase 1) possesses 
a body-centered tetragonal structure and is stable above 
414K. 1 The unit cell of this phase is depicted in Fig. T 
and consists of a lattice of SO4 tetrahedra, each bonded 
to a hydrogen via an O-H chemical bond. Each chem- 
ically bonded hydrogen also forms an O- • - H hydrogen 
bond with an oxygen of a neighboring SO4 tetrahedron. 
The resulting hydrogen-bond network becomes dynamic 
above the transition temperature and can visit a number 
of distinct topologies, owing to four possible oxygen bind- 
ing sites for each SO4 node in the network. The reigning 
view in the literatur o 9 ' 10 ' 11 ' 12 ' 13 is that long-range proton 
transport in superprotonic CSHSO4 occurs as the net re- 
sult of two separate mechanisms: first, the reorientation 
of the hydrogen-bond network by rapid rotations of the 
sulfate tetrahedra; and second, the hopping of the pro- 
ton between oxygens of neighboring tetrahedra across the 
O-H- • • O complex. The second step has generally been 
considered rate limiting and is thought to occur at fre- 
quencies of the order 10~ 9 s _1 , whereas the first is ex- 
pected to happen more frequently by at least two orders 
of magnitude. 14 ' 15 This two-step process is often referred 
to collectively as the Grotthus mechanism. 



II. METHOD 

We performed Car-Parrincllo molecular dynamics 
simulations^ of Phase-I CSHSO4 in the canonical NVT 
ensemble at temperatures of 550 K, 620 K, and 750 K (su- 
perheated), with temperatures maintained by means of 
a Nose-Hoover thermosta t 17 ! 18 ' 19 . Each simulation cov- 
ered 25 ps of thermalized dynamics following 5 ps of 




FIG. 1: (Color online) Structure of the conventional unit cell 
of Phase-I CsHSCk. Hydrogen atoms are shown in white, 
oxygen in red, sulfur in yellow, and cesium in blue. Hydrogen 
bonds are denoted by broken green lines. 



equilibration. This length of time proved sufficient for 
sampling several hundred to a thousand jump events, 
making stastical inferences possible. In each case, our 
supercell was comprised of 112 atoms (sixteen complete 
CSHSO4 units). Simulations were performed in a plane- 
wave basis using standard- valence ultrasoft pseudopoten- 
tials for hydrogen, oxygen, and sulfur; a 6s°' 5 5d 005 6p 005 
norm-conserving pseudopotcntial with nonlinear core 
correction for cesium; and the Perdew-Burke-Ernzcrhof 
exchange-correlation functional^. All pscudopotentials 
are obtainable from the Quantum-ESPRESSO website^. 
Cutoffs of 25 Ry and 150 Ry were used for the wavefunc- 
tions and charge density, respectively. The fictitious Car- 
Parrinello mass was /i = 700 with At = 7.5 au, and the 
lattice parameter was chosen based on the experimental 
value just above the superionic transition temperature, 
taken from Ref. I22I . 

In presenting the results of our simulations, we have 
divided the dynamics into categories of chemical-bond 
dynamics and hydrogen-bond dynamics. We include in 
our definition of chemical-bond dynamics any breaking or 
forming of O-H chemical bonds by Grotthus-type hop- 
ping of a proton between oxygens of neighboring tetra- 
hedra. Any change in the hydrogen-bond network struc- 
ture resulting from breaking or forming O-H hydro- 
gen bonds that does not also involve breaking or form- 
ing O-H chemical bonds is considered in the category of 
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FIG. 2: (Color online) Schematic depiction of (a) a sequence 
of chemical-bond jumps nucleated by the formation of an 
H2SO4 defect in the uppermost tetrahedron; (b) a hydrogen- 
bond network change induced by rotation of a host SO4 tetra- 
hedron; and (c) a hydrogen-bond network change resulting 
from a direct hydrogen-bond hop with little or no rotation of 
the host SO4 tetrahedron. The color scheme follows that of 
Fig. [T] with final configurations in jumping events shown as 
semi-transparent . 

hydrogen-bond dynamics. For additional clarity, Fig. [2] 
offers schematic representations of sample jump events 
from each category of dynamics. 



III. CHEMICAL-BOND DYNAMICS 

An O-H chemical bond was defined by considering 
interactions with oxygens within a cutoff distance of 
Ron < 1-15 A, which represents the initial separation be- 
tween the first and second coordination peaks of the cal- 
culated oxygen- hydrogen radial pair distribution function 
(RDF), displayed in Fig. [3] Classification as chemical or 
hydrogen bond proved more difficult for O-H pairs sepa- 
rated by an intermediate distance 1.15 < Ron < 1-35 A 
due to an inherent difficulty in resolving the overlap in 
the first two RDF peaks in that range. The ambigu- 
ity is also noticeable upon examination of the coordi- 
nation number ?i(r), which is nearly flat in this region. 
For such O-H pairs, we instead employed a history- 
dependent definition, basing the bond category on the 
last visited unambiguous bonding regime. Under this def- 
inition, an existent O-H chemical bond was considered 
broken only when Ron > 1-35 A; analogously, an exis- 
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FIG. 3: The oxygen-hydrogen radial pair distribution func- 
tion goti(r) (solid line, left axis) and the coordination number 
n(r) (dashed line, right axis), calculated from a simulation at 
620 K. Regions classified as chemical and hydrogen bonds are 
delineated, as well as the intermediate region for which the 
history-dependent definition was employed. 




FIG. 4: Autocorrelation functions for existence of (a) chemi- 
cal and (b) hydrogen bonds after a time t for simulations at 
550 K, 620 K, and 750 K, as well as at 620 K with fixed sulfate 
tetrahedra. 



tent O • • • H hydrogen bond was considered broken only 
when i? H < 1.15 A. 

Within the Grotthus mechanism, local proton trans- 
fer via a scries of correlated jumps prompts changes in 
the chemical-bond structure. Such jumps are first nu- 
cleated by the formation of a metastable H2SO4 defect, 
which subsequently propagates along the network back- 
bone, acting as a successive proton donor for neighboring 
tetrahedra at each stage. The individual jumps are them- 
selves short hops across a double-well potential barrier, 
where the two stable minima represent the O- • H and 
O-H distances and are separated by about 0.5 A. A sin- 
gle Grotthus hop therefore has the effect of swapping a 
chemical and a hydrogen bond, an action that is repeated 
as the proton transfer propagates across the hydrogen- 
bond network chain. This model of local proton transfer 
in superprotonic CSHSO4, represented schematically in 
Fig- 121 a), is easily observable in our simulations. Indeed, 
just over half (51%) of the chemical-bond jumps that we 
register at 620 K occur as a direct result of H2SO4 de- 
fect formation by the donation of a second proton from 
a neighboring tetrahedron, in accordance with the Grot- 
thus model. The remaining jumps are nucleated as a 
result of random local fluctuations in the bond structure. 

It is a straightforward process to track bond formation 
and annihilation, and we can define a time autocorrela- 
tion function for bond existence as C e (t) = (a(O)a(t)), 
where a(t) is 1 if a particular type of bond exists be- 
tween a given hydrogen-oxygen pair at time t, and 
otherwise. In practice, we can improve our statistical 
sampling by averaging C e (t) over all available time inter- 
vals of length t in the simulation. Using this quantity, 
we can obtain a detailed picture of the timescales of the 
chemical- and hydrogen-bond dynamics. The hydrogen- 
bond and chemical-bond existence autocorrelation func- 



tions are displayed in Fig. 0J 

Beyond about 20 fs, we observe a slow exponential de- 
cay in the chemical-bond existence autocorrelation C e (t), 
with characteristic 1/e decay times in the 11-15 ps range. 
These values are recorded in Table U along with average 
overall statistical frequencies of chemical-bond jumps for 
each temperature and the fraction of such jumps that 
subsequently reverse themselves. Here a chemical-bond 
jump is defined as a complete exchange of a chemical 
and a hydrogen bond across an O-H- • • O complex. The 
surprising commonality of chemical-bond jump events is 
reconcilable with the slow decay rate only when one con- 
siders the extremely high rate of jump reversal, which av- 
erages around 85% and features no significant variability 
with temperature. The combination of frequent jump- 
ing and a high, temperature-independent reversal rate 
is suggestive of a potential energy surface featuring an 
especially shallow activation barrier. 

At very short times (<10 fs), we observe a fast fall-off 
before transitioning to the slower decay regime. In this 
region, we are below the timescale of any jump rever- 
sal subsequent to chemical-bond breaking and forming, 
resulting in a much more rapid decay. The absence of 
any noticeable high-frequency periodicity in Fig. HJa) in- 
dicates that jump reversal carries no preferred timescale. 
Rather, it is likely that the reversal probability is a con- 
sequence of stabilization or destabilization of the local 
potential energy surface from SO4 tetrahedral reorienta- 
tions. 

As an indicator of the potential effect of the motion of 
the SO4 tetrahedra on the chemical- and hydrogen-bond 
dynamics, we have also chosen to run a second simulation 
at 620 K in which all ions except for the hydrogens were 
immobilized (denoted "fixed-S04" in Fig. HI) • The config- 
uration was chosen from a well-equilibrated timestep of 
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TABLE I: Various quantities derived from a statistical anal- 
ysis of the chemical-bond dynamics at 550 K, 620 K, and 
750 K. Featured columns indicate (1) the characteristic decay 
time r in an exponential fit Ae~ (t ^ T ^ of the long-time data 
in Fig. 21 a); (2) the average overall observed frequency v c of 
chemical-bond jump events per ion; and (3) the fraction of 
these events that subsequently reverse themselves. 



Temperature (K) 


T (ps) 


v c (THz) 


% rev. 


550 


15.2 


0.58 


83 


620 


12.5 


0.80 


85 


750 


11.2 


0.94 


85 



the fully mobile simulation. Interestingly, all chemical- 
bond dynamics ceased in this simulation. This fact is par- 
ticularly notable in light of previous investigation al 5 1 1 1 
which postulated that the most important factor in in- 
ducing chemical-bond jumping is a reduction in the 0-0 
distance across the O-H- • ■ O complex due to SO4 reori- 
entation. However, in our fixed-S04 simulation, not a 
single chemical-bond jump event was registered, despite 
the continuous presence of 0-0 distances as short as 
2.30 A across O-H- • • O complexes. Notably, this value is 
approximately 0.15 A shorter than the average 0-0 dis- 
tance across O-H- ■ • O complexes of hydrogens actively 
involved in chemical-bond jumps in the fully mobile simu- 
lation. Moreover, in our fully mobile simulations, we find 
only a 0.01 A difference between the average 0-0 dis- 
tances across O-H- • • O complexes of hydrogens involved 
in chemical-bond jumping events and ordinary hydrogens 
not involved in jump events of any sort. We therefore 
conclude that the primary contribution of the SO4 tetra- 
hedra to chemical-bond jumping results from their vibra- 
tional or rotational dynamics rather than simply their in- 
stantaneous orientation. Also, despite the low barrier for 
Grotthus-style chemical-bond hopping, dynamic degrees 
of freedom connected solely to the hydrogens are never- 
theless insufficient to permit chemical-bond breaking or 
forming, detailing the necessary role of the oxygen modes 
in that process. 

Our results reveal that the general chemical-bond jump 
frequency is relatively high; moreover, it is of the same 
order as the hydrogen-bond dynamics (compare Tables U 
and [IIJ . This contrasts with the view of chemical-bond 
jumping as substantially rate limiting and differing in 
timescale from the hydrogen-bond dynamics by two or 
more orders of magnitude. We instead find that the lim- 
iting factor in the effective rate of chemical-bond jumps 
is the extraordinarily high rate of jump reversal, which 
we suggest is linked to the dynamics of the SO4 tetrahe- 
dra. Yet even when jump reversals are considered, our 
effective chemical-bond dynamics are significantly faster 
than the proposed nanosecond scale. Instead, our results 
are consistent with recent experiment o 23 ' pointing to 
much faster chemical-bond dynamics than have hitherto 
been assumed, perhaps even on the picosecond scale. 

It should be noted that in our analysis of jump re- 



versals, we have considered only single jumps that are 
subsequently reversed. In so doing, we have neglected 
the reversal of collective sequences of jump events. Such 
higher-order jump sequences are especially difficult to 
properly consider given the inherent time- and length- 
scale limitations imposed by the first-principles method- 
ology. It is expected that accounting for the potential 
reversal of such longer sequences would limit the num- 
ber of counted "successful" jumps (those contributing to 
macroscopic proton diffusion). The magnitude of any 
such limitation is difficult to predict, however. 



IV. HYDROGEN-BOND DYNAMICS 

For the purposes of this work, we have defined a hydro- 
gen bond in terms of the oxygen- hydrogen distance alone, 
with the additional restriction that hydrogen bonding 
cannot involve oxygens attached to the same host sulfate 
tetrahedron as the hydrogen. Under this definition, the 
usual practice of restricting ZO-H- • • O had no apprecia- 
ble effect on the counted hydrogen bonds and was omit- 
ted for sake of simplicity. A hydrogen-bond maximum 
cutoff distance of Ron < 2.23 A was chosen based on the 
distance for which the coordination number n(r) = 2, in- 
dicating the tail end of the second coordation peak (as- 
sociated with H- • • O) in the calculated oxygen-hydrogen 
RDF (Fig. EJ). The minimum cutoff of i? H > 1.35 A 
was chosen based on the clear point of separation for the 
second RDF peak and the end of the plateau region in 
the coordination number. As for the chemical bonds, we 
implement a history-dependent definition for categorizing 
bonds in the intermediate range of 1.15 < Ron < 1.35 A. 

Fig. 0|b) shows the bond-existence autocorrelation 
function C e (t) for the hydrogen bonds in simulations at 
550 K, 620 K, 750 K, and for the "fixed-S0 4 " simulation 
at 620 K in which all ions except the hydrogens are im- 
mobilized. At longer times, we observe an exponential 
decay of the hydrogen bonds for the fully mobile simula- 
tions that far outpaces that of the chemical bonds. The 
graph also reveals that at short times (< 50 fs; see figure 
inset), the hydrogen bond network in the fixcd-S04 sim- 
ulation remains very dynamic, approximately following 
the equivalent curve for the fully mobile system. How- 
ever, C e (t) soon begins to oscillate around a fixed running 
average, indicating repeated visitation of a few alternat- 
ing configurations. Interestingly, the overall frequency of 
hydrogen-bond breaking is actually greater for the fixed- 
SO4 simulation than for the fully mobile simulation. 

Table [TT] contains these hydrogen-bond breaking fre- 
quencies, as well as likelihoods for reversal of hydrogen- 
bond network reorganization events. In addition to pro- 
viding overall values, we have divided the hydrogen-bond 
dynamics into two categories based on the location of 
the newly formed hydrogen bond with respect to its pre- 
decessor. Our first category consists of hydrogen bonds 
transferred between oxygens of the same destination SO4 
tetrahedron; the second contains hydrogen bonds trans- 
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ferred between oxygens of neighboring SO4 tetrahedra. 
In practice, higher temperatures generally show a greater 
preference for exchanges between oxygens of different 
tetrahedra than do lower temperatures, but in all cases, 
such exchanges outnumber those between oxygens of the 
same tetrahedron by a margin of two or three to one. 
Fixing the SO4 tetrahedra pushes that margin even fur- 
ther. 

Unlike in the case of the chemical-bond dynamics, 
freezing the degrees of freedom of the SO4 tetrahedra 
does not prevent reconfiguration of the hydrogen-bond 
network via bond breaking and forming. In fact, Ta- 
ble |TT] reveals that inhibiting SO4 rotation actually en- 
hances the frequency of hydrogen-bond breaking and 
forming, particularly for bonds exchanged between oxy- 
gens of neighboring tetrahedra. However, some degree of 
rotation is required in order to visit a larger region of con- 
figuration space and prevent repeated visitation of iden- 
tical configurations, a vital stipulation for macroscopic 
proton transport. 

The point of separation of the hydrogen-bond auto- 
correlation curves for the fixed-S04 and fully mobile sys- 
tem, as seen in the inset of Fig. [31(b) , can be interpreted 
physically as a characteristic reversal time: if a newly 
formed hydrogen bond is to be accepted, enough SO4 
rotation must occur within the approximately 50 fs win- 
dow to sufficiently alter and imbalance the energy land- 
scape, thereby minimizing likelihood of back hopping. 
Accordingly, Table [TT] shows that the observed fraction of 
hydrogen-bond network reorganization events that sub- 
sequently reverse themselves within this time period in 
the fixed-S04 simulation is more than double that of the 
fully mobile simulation at the same temperature (38% 
versus 77%). 

A Fourier transform of the hydrogen-bond existence 
autocorrelation function C e (t) gives a good measure of 
the typical oscillation frequencies for the hydrogen-bond 
forming and breaking in the fixed-S04 simulation. Fig. [5] 
compares this result to the vibrational density of states 
for the hydrogens in that simulation as well as in the fully 
mobile simulation. A comparison of Figs. E^b) and (c) 
allows us to distinguish the hydrogen vibrational modes 
that are not directly connected to hydrogen-bond break- 
ing from those that are. Those not linked to changes 
in the hydrogen-bond network are represented by clus- 
ters of broader peaks around 15-20 THz, 30-40 THz, 
and 80-95 THz. Of these, the two lowest-frequency clus- 
ters most likely represent bending modes, whereas the 
highest-frequency cluster contains stretching modes. The 
primary peaks associated with bond breaking and form- 
ing are a low-freqency peak near 9 THz and a fundamen- 
tal second peak at 28 THz, along with its accompanying 
overtone peaks at higher frequencies. Since these peaks 
are completely suppressed in the result for the fully mo- 
bile simulation shown in Fig. [5^a), they represent the 
signature oscillations inhibited upon stabilization of new 
configurations by reorientations of the SO4 tetrahedra. 
Notably, the half-period switching time represented by 



(a) H vibrational 




H vibrational (fixed S0 4 ) 


(c) 


HB breaking (fixed S0 4 ) 
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Frequency (THz) 

FIG. 5: Vibrational density of states for the protons in (a) 
the fully mobile and (b) the fixed-SCU simulations, along with 
(c) the hydrogen-bond switching frequency spectrum for the 
fixed-S04 simulation. Data are from simulations at 620 K, 
and densities of states were obtained from a Fourier transform 
of the appropriate velocity autocorrelation function. 

the low-frequency peak matches the characteristic rever- 
sal threshold obtained from Fig. [4] The higher-frequency 
peak at 28 THz is also evident in that same figure, ap- 
pearing as shallow oscillations at short timescales. The 
existence of pronounced overtones for the 28 THz peak in 
both the vibrational and bond-breaking frequency spec- 
tra suggests significant anharmonicity in the potential for 
the H- • ■ O bond. 

Since the primary effect of suppressing SO4 rotation is 
to encourage reversal of hydrogen-bond network reorien- 
tation phenomena rather than to inhibit such reorienta- 
tions altogether, SO4 rotation alone cannot satisfactorily 
account for the full network dynamics. Rather, we ob- 
serve that the dominant mechanism for hydrogen-bond 
network reorganization involves hydrogen-bond transi- 
tions that are best described as rapid, discrete angu- 
lar jumps between two stable states rather than as a 
smooth evolution driven by SO4 tetrahedral rotations, as 
has generally been proposed previously. These two states 
correspond to different orientations for which the ZS-O- 
H angle for chemically bonded hydrogens is maintained 
near the tetrahedral angle of 109.5°. This phenomenon 
is depicted schematically in Fig. [5Jc) and resembles the 
proposed diffusion mechanism in a recent simulation of 
liquid water—. 

Additional evidence for the angular hopping model ap- 
pears in Fig. [6l which outlines distributions for certain 
geometrically relevant angles in both the fully mobile and 
fixed-S04 simulations at 620 K. The ZS-O-H angles for 
chemically bonded protons have a relatively small spread 
and are peaked around the described tetrahedral geom- 
etry. Angles greater than 145° are not represented, in- 
dicating that the protons lie primarily on the surface of 
a cone centered on the S-0 bond and with a half-angle 
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TABLE II: Various quantities derived from a statistical analysis of the hydrogen-bond dynamics at 550 K, 620 K, 750 K, and 
for the fixed-SC>4 simulation at 620 K. The data is divided into statistics for hydrogen-bond exchanges between oxygens of 
the same SO4 tetrahedron, those of different SO4 tetrahedra, and overall totals for either type of exchange. Featured columns 
include (1) the fraction of total hydrogen-bond exchanges representing a particular class of exchange; (2) the average overall 
observed frequency vu of hydrogen-bond jump events per ion; and (3) the fraction of these events that reverse themselves within 
a 50 fs window. 







Same SO4 






Different S0 4 




Overall 




Temperature (K) 


% tot. 


Vk (THz) 


% rev. 


% tot. 


»b (THz) 


% rev. 


u b (THz) % rev. 


550 


38 


0.47 


49 


62 


0.79 


34 


1.16 


40 


620 


25 


0.58 


38 


75 


1.72 


38 


2.30 


38 


750 


23 


1.01 


35 


77 


3.45 


39 


4.46 


39 


620 (fixed-S0 4 ) 


21 


0.96 


67 


79 


3.57 


79 


4.53 


77 




45 90 135 180 45 90 135 180 



ZS-O-H (deg) 

FIG. 6: S-O-H angles for (a) chemically bonded and (b) 
hydrogen-bonded protons in fully mobile (solid) and fixed- 
SO4 (dashed) simulations at 620 K. 

of 65-70°. Notably, the angular distribution does not 
change appreciably between the fully mobile and fixed- 
SO4 simulations. Although hydrogen-bonded protons are 
generally less constrained, the ZS-0--H distributions 
for both simulations are still peaked near 110°. However, 
for the fixed-S04 simulation, a second peak appears in 
the angular distribution at around 140° as a byproduct 
of the hydrogen-bond hops. The separation of the two 
peaks for the fixed-S04 case in Fig. [6^b) indicates that 
generally a net SO4 tetrahedral reorientation of around 
30°, involving either the hydrogen-bond donor or accep- 
tor tetrahedron, accompanies the hop to alleviate the 
lattice strain it induces. This value for the SO4 angu- 
lar rotation agrees well with what has been proposed in 
the literatur o 10 i 22 ' 26 i 27 , but our resolution is insufficient 
to pinpoint which of the particular competing models 
is most likely to be correct. Such reorientation is also 
responsible for altering the potential energy surface to 
minimize back hopping. 



V. ROTATION OF SO4 TETRAHEDRA 

Since it has already been established that the 30° reori- 
entation of the SO4 tetrahedra must take place within a 
50 fs window to maximize the potential for non-reversing 
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co (rad/ps) 

FIG. 7: (a) Angular velocity profile for the SO4 tetrahedra at 
620 K (solid) , along with the corresponding fraction of tetra- 
hedra n(ui) with velocities < ui (dashed). Inset (b) shows the 
angular distribution of the sulfate tetrahedral axes of rotation 
u! with respect to tsh, and inset (c) shows the radial distri- 
bution of the u> vectors in space with respect to the primary 
crystallographic axes. 



transitions, we can obtain 10.5 rad/ps as a back-of-the- 
envelope estimate for the minimum SO4 angular velocity 
required to prevent reversal following a hydrogen-bond 
switch. The SO4 angular velocities follow a Boltzmann 
distribution and are plotted in Fig. [7] The plot reveals 
that velocities of this magnitude, although somewhat 
rare, are indeed accessible, representing about 8-9% of 
the SO4 tetrahedra in the 620 K simulation at any given 
time. 

Fig- IZIc) shows an isosurface of the angular velocity 
unit vectors Cj for the rotation of the SO4 tetrahedra, 
averaged over all such groups in the 620 K simulation. 
Areas of high density therefore represent preferred axes 
of rotation, a clear structure for which is visible in the fig- 
ure. These rotation axes do not align towards the chem- 
ically bonded hydrogen or its accompanying oxygen, as 
is evident from Fig. [7](b) . Instead, they orient along the 
edges of a cube rotated 7r/4 in the (001) crystal plane 
with respect to the conventional unit cell, thereby cor- 
relating with the centers of the nearest-neighbor tetra- 
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FIG. 8: Distribution of angular distances traveled by sulfate 
tetrahedra as a function of time at 620 K. Successive curves 
represent values at 0.5, 1, 2, 3, and 5 ps. The inset gives a 
measure of the fraction of tetrahedra with greater than 70° 
rotation from their initial positions as a function of time at 
550 K, 620 K, and 750 K. 



hedra. The SO4 rotational orientations thus appear to 
be governed by the locations of nearby SO4 tetrahedra 
rather than the location of the locally bonded hydrogen 
or the orientation of its corresponding chemical bond. In 
other terms, rotational axes align along ?ss rather than 
? SH or f so- 
In view of the fact that angular jumps alone, stabi- 
lized by small, rapid tetrahedral reorientation events of 
about 30°, account for much of the hydrogen bond net- 
work dynamics, it is desirable to analyze and requantify 
the relative contribution of larger-scale rotational mo- 
tion. Fig. [3 which shows the averaged distributions of 
angular distances traveled by SO4 tetrahedra as a func- 
tion of time, illustrates one effect of slower rotational dy- 
namics. There is evidence of the appearance of a second 
peak in the distributions representing a new, stable equi- 
librium configuration at 75-80° rotation with respect to 
the original tetrahedral orientation. This peak begins to 
manifest in statistically measurable quantities only after 
250 fs (see inset of figure) , making the quickest of such re- 
orientation events several times slower than the timescale 
of the fast 30° reorientation event described previously. 
We note that the faster dynamics with the smaller reori- 
entation angle cannot be distinguished in Fig. [5] since it 
is buried within the first peak, which primarily depicts 
the fast librational modes. 

Additional confirmation of the existence of multiple 
distinct timescales for the tetrahedral orientation can be 
seen in Fig. [9l which portrays the angular time autocor- 
relation of tetrahedral configurations. We calculate this 
quantity according to Cg (t) = (fso(i)'fso(O)), producing 
a measure of the cosine of the average angular distance 
traveled by an S-0 unit vector f so as a function of time. 
As before, we maximize our available statistics by aver- 
aging Cg(t) over all available time intervals of length t 
in the simulation. The slower rotation events manifest 



FIG. 9: Autocorrelation function for angular distance traveled 
between SO4 tetrahedral orientations separated by time t. 



themselves as a quasi-linear decay in the autocorrelation 
at longer times. At short times (< 250 fs), the slope is ap- 
preciably steeper, indicating faster dynamics on average. 
Expectedly, the separation between these two regimes 
agrees with the timescale of the emergence of the slow 
rotation in Fig. [5] In addition, a shoulder indicating the 
timescale of a rotation to a nearby local minima is clearly 
distinguishable at around 50-60 fs, in agreement with our 
previous indicators of fast dynamics on that scale. This 
shoulder repeats itself as periodic oscillations and is also 
detectable as a peak in the Fourier transform (not shown) 
of the curve in Fig. [Sj The oscillations also span the in- 
termediate region in which rotations of both short and 
long timescales are manifest before they are lost in the 
slower dynamics at longer times. We further note that 
changes in temperature have no appreciable effect on the 
timescale of the fast rotation, as measured by the loca- 
tions of the shoulders and oscillations in the figure. This 
is consistent with a picture in which the reorientation is 
connected to a hydrogen-bond hop and therefore has an 
almost negligible rotational barrier. 

The anisotropy of the diffusion tensor for CSHSO4 
has been documented experimentally^ and is a geo- 
metric consequence of the greater angular distance that 
must be traversed by a diffusing proton traveling across 
two sulfate layers along the [001] or [00T] directions 
(117°) compared to a similar journey in a {001} plane 
(78°). Our results confirm that diffusion parallel to the 
(001) plane dominates: the directional mean-square dis- 
placement of the hydrogen atoms rises 2.5 to 5 times 
faster along the [010] or [100] directions than along the 
[001] direction, with higher temperatures favoring higher 
anisotropy. This disparity suggests a different dominant 
mechanism for diffusion along a (001) direction, in part 
because the angular distance is too great to be easily 
accommodated by the described hydrogen bond hopping 
mechanism, and in part because the corresponding angu- 
lar velocities that would have to be reached by the sulfate 
tetrahedra to prevent backhopping in such a scheme are 
unreasonably high. Instead, we find that slower SO4 ro- 
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tation plays the dominant role in overcoming the larger 
barrier, in line with the more traditional model of proton 
transport in CSHSO4 (depicted in Fig.^b)). In this case, 
the dynamics are slow enough that the timescale of the 
large-scale rotation is no longer a hindrance. We further 
note that the approximately 40° difference in the angular 
distance that must be traveled in the [001] direction with 
respect to a similar journey in the (001) plane, combined 
with the aforementioned 30° SO4 reorientation, satisfac- 
torily accounts for the appearance of the second peak in 
Fig- El to within a rough estimate. This suggests that for 
diffusion in the [001] direction, the slow rotation is prob- 
ably followed by a rapid hydrogen-bond hop, although 
the former clearly determines the timescale. 



VI. HYDROGEN-BOND NETWORK 
TOPOLOGY 

We have also analyzed the basic topology of the 
hydrogen-bond network. Table IIIII shows the relative 
probabilities of various SO4 bonding configurations at 
550 K, 620 K, and 750 K, organized according to the 
number of hydrogen bonds donated (Nd) and accepted 
(N a ) by the SO4 tetrahedron. In our definition, a 
donated bond is formed between a local chemically 
bonded hydrogen and an oxygen on a neighboring tetra- 
hedron; an accepted bond is a hydrogen bond formed 
between a local oxygen and a hydrogen that is chemi- 
cally bonded to a neighboring tetrahedron. For an ideal 
one-dimensional network, one would expect all tetrahe- 
dra to have Nd=N a =l. In our simulations, we observe 
this type of ordinary link only 41-61% of the time, sug- 
gesting the actual network topology is much more com- 
plicated. Other highly probable configurations include 
one with (N a , Nd) — (1, 0), which can be thought of as a 
terminator in the hydrogen bond network (18-26%); and 
one with (N a ,Nd) = (1,2), which can be interpreted as 
a network branching point (8-14%). According to the 
table, there is a great deal of variability in the number of 
hydrogen-bond donors Nd, but configurations with N a j^l 
are comparatively rare. The primary effect of increasing 
the temperature seems to be a decrease in the number of 
ordinary linear network links with Nd=N a =l in favor of 
network terminators with (N a ,Nd) — (1,0), resulting in 
a more nodal network. Some of the rare (but nonetheless 
statistically significant) configurations are signatures of 
Grotthus-type jumps in progress: immediately following 
a standard chemical-bond jump from one SO4 tetrahe- 
dron to another, nucleated at a link in an ordinary lin- 
ear chain, the source tetrahedron registers a topological 
configuration of the form (N a , Nd) = (2,0), whereas the 
H2SO4 destination tetrahedron acquires a configuration 
of the form (N a ,N d ) = (0,2). 

Table ITVl lists the average lifetimes of the local topolo- 
gies listed in Table IIIII Although these values are av- 
erages and do not account for the complete distribution 
of possible lifetimes as do the autocorrelation curves of 




FIG. 10: Existence autocorrelation functions for four of 
the most likely tetrahedral bonding configurations. In the 
schematic diagrams, a solid line represents a chemical bond, 
a dashed line represents a hydrogen bond, and an "X" indi- 
cates the absence of a bond. 



Fig. [TU1 they are nonetheless useful for purposes of quali- 
tative comparison. Lifetimes are generally well correlated 
with relative frequencies, with hydrogen bonding to a sin- 
gle secondary tetrahedron (N a =l) acting as a stabilizing 
force. As the temperature increases, the lifetimes of con- 
figurations with N a ^l are affected very little, but we 
observe a sharp systematic decline in nearly all configu- 
rations with N a =l. Notably, this trend does not always 
follow that of the relative frequencies in Table IIIII For 
example, network terminators with (N a ,Nd) = (1,0) ex- 
hibit a decrease in average lifetime but an increase in 
their overall commonality, suggesting a more nodal but 
also more dynamic network at high temperatures. 

Tables IIIII and IIVI neglect hydrogen bonds to the same 
SO4 tetrahedron as the host, in accordance with our orig- 
inal definition of the hydrogen bond. However, if we re- 
lax the hydrogen-bonding restriction requiring bonds to 
be between oxygens of different tctrahcdra, we find that 
these "self-hydrogen bonded defects" , although short- 
lived, are nonetheless relatively common, representing 
about 3% of the total O- • • H interactions at 620 K. More- 
over, we find that hydrogens involved in chemical and hy- 
drogen bonds to the same SO4 tetrahedron do not permit 
simultaneous hydrogen bonding to oxygens of neighbor- 
ing tetrahedra, meaning these complexes function as ter- 
minators for the hydrogen-bond network chains. Also, 
the average velocity of oxygens in self-hydrogen bonded 
SO4 complexes is consistently about 10% higher on av- 
erage than in ordinary complexes with Nd=N a =l. This 
suggests that underbonding lessens the degree of con- 
straint and enhances oxygen mobility in such units, likely 
aiding further reorganization of the hydrogen-bond net- 
work. 

In addition to examining the network topology in 
terms of connectivity between neighboring SO4 tetrahe- 
dra, one can obtain a slightly different topological gauge 
by looking at the number of chemical and hydrogen bonds 
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TABLE III: Observed relative probabilities (%) of SO4 bonding configurations, organized according to the number of hydrogen 
bonds accepted (rows) and the number donated (columns) by the SO4 tetrahedron. Data are from simulations at 550/620/750 K. 



% Prob. 


N d = 


N d = l 


N d = 2 


N d = 3 


N a = 


2/3/7 


3/4/6 


1/1/2 


1/1/1 


N a = 1 


18 / 22 / 26 


61 / 53 / 41 


11 / 15 / 14 


1/1/1 


N a = 2 


1/1/2 


<1/<1/1 


0/0/0 


0/0/0 



TABLE IV: Average lifetimes (fs) of SO4 bonding configurations, organized according to the number of hydrogen bonds accepted 
(rows) and the number donated (columns) by the tetrahedron. Data are from simulations at 550/620/750 K. 



Lifetime 


(fi>) 


N d = 


N d = 1 


N d = 2 


N d = 3 


N a = 





48 / 53 / 55 


47 / 43 / 42 


21 / 19 / 25 


34 / 21 / 24 


N a = 


1 


176 / 159 / 117 


238 / 169 / 105 


128 / 119 / 75 


81 / 92 / 54 


N a = 


2 


27 / 20 / 24 


20 / 17 / 20 


0/0/0 


0/0/0 



formed by a single proton. Fig. [TO] shows the existence 
autocorrelation curves for four of the most common topo- 
logically distinct bonding configurations for a proton. 
These curves give an idea of the characteristic decay 
times for each configuration. The result for a standard 
bonding configuration, in which a proton forms one chem- 
ical and one hydrogen bond (^4) is shown for timescale 
reference. We note that topologies with protons forming 
a single chemical bond but no hydrogen bond (B) are 
surprisingly stable, with a characteristic decay time that 
is relatively large on the timescale of both the hydrogen- 
bond hopping and its corresponding strain-relaxing SO4 
reorientation. Defects of this class are less constrained 
and more mobile than their ordinary counterparts, and 
as in the case of the self-hydrogen bonded complexes, 
the increased mobility facilitates network reorganization 
much more readily. Network-branching configurations 
with multiple hydrogen bonds (C) play a more direct 
role in the reorganization of the hydrogen-bond network 
but have only intermediate stability. Configurations with 
no chemical bond (D) are extremely short-lived. 

Techniques involving the graph-theoretic adjacency 
matrix (see Appendix) also offer a convenient way of 
characterizing the topology of the overall network and 
extracting configurations most likely to induce a diffu- 
sive event. In particular, we are able to further classify 
the network topology in terms of rings, meaning some 
part of the network ultimately connects back to itself 
in a closed loop; and chains, meaning the network ei- 
ther remains linear or branches, with the restriction that 
any two network vertices are connected by exactly one 
unique directed path (i.e., a graph-theoretic tree). The 
specifics of our classification algorithm are described in 
the Appendix. It should be noted that such a dichotomy 
requires classification of every node as either a chain or 
a ring but does not allow any given node to be doubly 
counted as belonging to both categories. Table [V] lists 
the likelihood of finding a tetrahedron in various ring 
and chain topologies in an ordinary simulation timestep 
versus a timestep immediately preceding a chemical- or 



hydrogen-bond jump event. 

For ordinary configurations not involving a jump event, 
the network favors rings over chains at 550 K, whereas 
the trend is reversed at 750 K. The intermediate temper- 
ature of 620 K is a topological transition zone and shows 
a marked increase in configurations simultaneously con- 
taining both rings and chains. There is also a notable 
decrease in the average length of a chain and a slight de- 
crease in the average size of a ring at 620 K compared 
to the other temperatures. This is a further indication 
that the network is midway in a transition process from 
primarily rings to primarily chains, as a network configu- 
ration containing both would tend to inhibit the growth 
of either one at the expense of the other. 

Our findings indicate that at all temperatures, the 
presence of topological chains has a dramatic effect in en- 
hancing the likelihood of occurence of either sort of jump 
event. Conversely, the presence of rings strongly inhibits 
jumping. The trend is much more pronounced when one 
examines frames containing either rings or chains as the 
only topological species: in frames preceding a chemical- 
bond jump, we see a 23-24% increase in likelihood of the 
frame to contain exclusively chains and a corresponding 
13-23% decrease in its likelihood to contain exclusively 
rings, as compared to a frame in an ordinary non-jumping 
configuration. This difference is especially pronounced 
at higher temperatures. Nonuniform configurations con- 
taining both rings and chains also show an overall de- 
crease in jump likelihood. The trends in the chain and 
ring data also evidence that from a purely topological 
perspective, a configuration favorable for a hydrogen- 
bond jump lies midway between an ordinary configura- 
tion and one favorable for a chemical-bond jump. We 
conclude that the network ring and chain topology is 
a good indicator of both hydrogen- and chemical-bond 
jump likelihood and is substantially more effective in that 
capacity than the measure of the oxygen-oxygen distance 
across the O-H- • • O complex. 

It is worthwhile mentioning that periodic boundary 
conditions and limited supercell sizes have two major 
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topological consequences that must be considered in any 
analysis: first, they decrease the maximum length of 
chains that can be formed; and second, they tend to arti- 
ficially inflate the number of smaller rings, since creation 
of periodic images tends to wrap the network back onto 
itself prematurely. As such, the values in Table W\ should 
not be taken as absolutes, but qualitative comparisons 
are nonetheless useful and relevant. 



VII. PROTON KINETICS AND THE ISOTOPE 
EFFECT 

We can estimate the general three-dimensional pro- 
ton self-diffusion coefficient from the mean-square dis- 
placement (MSD) using the Einstein relation: Dh = 
lim^oo gj (MSD(t)}. Using this method, we estimate Dh 
in our simulations to be 1.8-3.5xl0 -6 cm 2 /s over the 
550-750 K temperature range, about an order of mag- 
nitude greater than extrapolations of experimental mea- 
surements to the same range (1.9-3.4x 10 -7 cm 2 /s) 15 . 
However, despite the error in the magnitude of the cal- 
culated diffusion coefficients, we nonetheless observe the 
proper scaling of Dh with temperature, indicating cor- 
rect calculation of the energetic barriers for diffusion. 

Multiple factors may contribute to the discrepancy 
of the calculated diffusion coefficients with experiment. 
First, calculations of these coefficients are notoriously 
difficult to converge, particularly when diffusivities are 
small and statistics are limited. One possible gauge of 
accuracy in the calculation can be achieved by comparing 
the result to the same quantity found using the Green- 
Kubo relation. The Green-Kubo method is based on the 
velocity autocorrelation function rather than the mean- 
square displacement; however, this method generally ex- 
hibits even poorer convergence in the absence of excel- 
lent statistics. Instead, we estimate the accuracy of our 
diffusion coefficients based on twice the standard error 
in the determination of the mean-square displacement 
slope, which results in an error estimate of ±15%. 

Second, we have already mentioned the difficulty in 
accounting for correlation and reversal of collective se- 
quences of jump events. In the limit of the short length- 
and timescales accessible by first-principles methods, it 
is expected that any diffusive correlations that persist 
over large regions of either time or space would be lost, 
resulting in inflated diffusion statistics. Similarly, jump 
events that remain localized due to high incidences of 
back-hopping and correlation are counted as contribu- 
tors to diffusion on the timescale accessible to our simu- 
lations, whereas these would not appear in a macroscopic 
measurement. 

We have also mentioned the artificially high number 
of small rings due to periodic boundary conditions as a 
potential consequence of small supercell size. Table [V] 
suggests a possible correlation between smaller average 
ring size and a higher propensity to jump, particularly 
at the lower temperature, meaning our finite size effects 



could result in a measurable increase in jump statistics 
and diffusion coefficients. However, any such decrease in 
average ring size from the periodic boundary conditions 
is likely to be accompanied by an increase in the overall 
number of rings. Since we have established the overall 
propensity for ring existence to inhibit hydrogen- and 
chemical-bond dynamics, the competition between these 
two effects should attenuate any potential impact on the 
macroscopic properties. 

It should be noted that we are neglecting any quantum 
behavior of the protons. However, an analysis of experi- 
ments on CSDSO4 suggests the isotope effect is relatively 
smalP-i 2 ^. In addition, theoretical work on the topologi- 
cally similar material KDP 30 concluded that the predom- 
inant effect of quantum derealization of the protons was 
limited to structural considerations, in that it decreased 
the H- • • O-H distance and consequently also the lattice 
parameter. The KDP analysis is also consistent with ex- 
perimental comparisons of CSDSO4 and CSHSO 4 10 ' 22 . 

A closer examination of the specific rate-limiting mech- 
anisms covered in our analysis of proton diffusion in 
CSHSO4 provides additional insight into the lack of any 
substantial isotope effect. Our findings indicate that 
the chemical-bond dynamics are much faster than pre- 
vious analyses have suggested, of the same order as the 
hydrogen-bond dynamics. Rather, the primary limita- 
tion is manifest in the dynamics of the SO4 tetrahc- 
dra, since these appear to govern the reversal rates of 
both chemical-bond jumping and hydrogen-bond hop- 
ping. Accounting for proton tunneling across the O- 
H- • • O double-well potential would therefore have little 
effect on the overall jump statistics, since the mobility 
of the SO4 groups is classically controlled. In addition, 
we have established the importance of the chain and ring 
topology of the hydrogen-bond network in promoting dif- 
fusive events. However, changes in the topology are pro- 
moted by two factors — hydrogen-bond hopping and slow 
rotation of the SO4 tetrahedra. Slow SO4 rotation is 
clearly a classical phenomenon, and hydrogen-bond hop- 
ping is coupled to a 30° reorientation of the heavy SO4 
tetrahedron, meaning its dynamics are also ultimately 
classical. 



VIII. CONCLUSIONS 

We have presented a detailed analysis of proton dy- 
namics in superprotonic Phase-I CSHSO4 based on first- 
principles molecular dynamics simulations. Our results 
confirm that the chemical-bond dynamics are dominated 
by local Grotthus-style hops which propagate succes- 
sively along the hydrogen-bond network backbone. Indi- 
vidually, these hops are comparatively frequent, pointing 
to a low diffusion barrier, but the net effective rate of the 
chemical-bond dynamics is limited by an anomalously 
high rate of jump reversal. We find that the propensity 
for such forward- and back-hopping along the O-H- • • O 
complex is in turn heavily influenced by the dynamics of 



11 



TABLE V: Observed relative probabilities (%) for various hydrogen-bond network topologies in an ordinary timestep, compared 
with similar quantities for timesteps immediately preceding a chemical- or hydrogen-bond jump event. Also listed are the 
relevant average ring and chain sizes for frames where those topologies exist. Ring sizes are calculated based on the number of 
tetrahedra involved in the ring, whereas chain sizes denote the maximum individual branch length within the graph-theoretic 
tree structure. Data are from simulations at 550/620/750 K. 



Description 


No jump 


Hydrogen-bond jump 


Chemical-bond jump 


Contains rings 


71 / 78 / 52 


64 / 72 / 49 


47 / 55 / 28 


Contains chains 


60 / 74 / 70 


69 / 81 / 75 


83 / 87 / 84 


Contains only rings 


40 / 26 / 30 


31 / 19 / 25 


17 / 13 / 16 


Contains only chains 


29 / 22 / 48 


36 / 28 / 51 


53 / 45 / 72 


Contains both rings and chains 


31 / 52 / 22 


34 / 52 / 24 


29 / 42 / 12 


Average ring size 


5.6 / 4.5 / 4.7 


5.3 / 4.5 / 4.6 


4.9 / 4.5 / 4.5 


Average chain size 


6.9 / 5.4 / 6.2 


6.9 / 5.4 / 6.0 


6.8 / 5.6 / 6.0 



the SO4 tetrahedra rather than by static local geometry 
alone. 

We have also shown that the dynamics of the hydrogen- 
bond network are dominated by fast, discrete angular 
jumps between neighboring oxygens rather than by slow 
rotations of the SO4 tetrahedra. Such jumps occur with 
greater frequency between oxygens belonging to different 
SO4 tetrahedra than between oxygens of the same tetra- 
hedron, by a factor of two or three. The hydrogen-bond 
jumps are accompanied by an approximately 30° reorien- 
tation of the participating SO4 tetrahedra to alleviate the 
lattice strain induced by the hop, thereby minimizing the 
likelihood of jump reversal. We have isolated a window 
of 50 fs for successful completion of this "fast" reorienta- 
tion event and showed that it exists independently of a 
second, slower reorientation mechanism, operating on a 
timescale at least five times greater than its counterpart. 
The slower mechanism amounts to ordinary SO4 rotation 
on a longer timescale, and we propose this to be the dom- 
inant hydrogen-bond network reorientation mechanism 
for diffusion along the [001] direction, for which angu- 
lar hops are significantly more difficult and less frequent, 
owing to the anisotropy of the CSHSO4 lattice. 

Our topological analysis of the hydrogen-bond net- 
work revealed a significant number of branching and 
network-terminating nodes, indicating a substantial devi- 
ation from linearity, particularly at higher temperatures. 
We postulate that the underbound network terminators 
play a role in network reconfiguration by aiding SO4 ro- 
tational mobility. Graph-theoretic methodology offered 
a way to isolate chains and rings as dominant topological 
features in the network, and we discovered that the pres- 
ence of chains and the absence of rings is a substantial 
predictor of likelihood for either a hydrogen- or chemical- 
bond jump event to occur. We propose that our topo- 
logical analysis could be easily extended to similar well- 
defined, hydrogen-bonded network solids. 

Finally, we apply our analysis to offer an explana- 
tion for the lack of a significant isotope effect in the 
CSHSO4/CSDSO4 system. In particular, we tie both the 
chemical- and hydrogen-bond dynamics to the classical 



dynamics of the SO4 tetrahedra and argue that the in- 
clusion of proton quantum tunneling should play a rela- 
tively minor role in the rate-limiting steps of the diffusion 
mechanism. 
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X. APPENDIX: GRAPH-THEORETIC 
METHODOLOGY 

Here we offer a detailed account of the graph-theoretic 
methods used to calculate the hydrogen-bond network 
topology. We represent the network as a directed graph 
with the edge vector pointing along the 0-H- ■ • O bond 
direction; that is, from the sulfate tetrahedron acting as 
the hydrogen-bond donor in the complex to the sulfate 
tetrahedron acting as the hydrogen-bond acceptor. The 
adjacency matrix Aij can then be constructed as an N x 
N matrix, where N is the number of tetrahedra in the 
unit cell and the indices i and j run over the donor and 
acceptor sulfate groups, respectively: 

I 1 if there exists a direct link i —> j 
[0 otherwise 

The diagonal elements An are set to zero. Topological 
characterization of a single node is then a straightforward 
process of performing a row sum to get the number of 
nodes to which it donates (Nd) and a column sum to get 
the number of nodes donating to it (N a ). It is also easy 
to categorize jump events by analyzing the difference of 
the adjacency matrices of successive timesteps. 

To determine ring connectivity and size, we exploit the 
property of adjacency matrices^, that element of A n 
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gives the number of unique directed pathways from i to 
j of length n. We take the size of the ring containing 
the i th node to be the lowest value of n in the interval [2, 
(N-l)} for which the diagonal element ^ 0. If A% = 
for all n in the interval, the node is not considered part 
of a ring. 

Deriving chain sizes is somewhat more complex, since 
we must account for multiple branching topologies and 
for topological mixtures of rings and chains. We first 
decompose the network into clusters of unconnected sub- 
graphs. This is done using the connectivity matrix Cy, 
which has the property that dj—1 if there is a path of 
any length connecting nodes i and j. We form the sym- 
metric Cij from Aij using WarshalPs algorithm 32 . 

We proceed to construct a matrix Sij that contains the 
shortest path between each pair of nodes i ^ j, by 

taking the lowest value of n in the interval [2, (A^-l)] for 



which A?, ^ 0. We can then take the maximum value of 
Sij across the columns to get a row array of maximum 
path lengths for chains originating at the i th node. The 
chain size is then determined by finding the maximum 
value of the resultant array over the nodes contained in 
each connected cluster, which can be easily determined 
by parsing dj . Finally, we add the restriction that none 
of the links in the chain can themselves be members of 
rings. This prevents counting of chains that are ficti- 
tiously long due to intermediate or terminating rings and 
ensures a clear separation between ring and chain topolo- 
gies. The resulting chain size is then calculated as the 
maximum span of the graph-theoretic tree. 

Very large networks may necessitate more efficient al- 
gorithms due to the expense of calculating A N ^ . How- 
ever, the system sizes in our study are sufficiently small 
to readily allow calculations using the described method. 
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